Matter-wave dark solitons: stochastic vs. analytical results 
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The dynamics of dark matter-wave solitons in elongated atomic condensates are discussed at finite 
temperatures. Simulations with the stochastic Gross-Pitaevskii equation reveal a noticeable, exper- 
imentally observable spread in individual soliton trajectories, attributed to inherent fluctuations 
in both phase and density of the underlying medium. Averaging over a number of such trajec- 
tories (as done in experiments) washes out such background fluctuations, revealing a well-defined 
temperature-dependent temporal growth in the oscillation amplitude. The average soliton dynamics 
is well captured by the simpler dissipative Gross-Pitaevskii equation, both numerically and via an 
analytically-derived equation for the soliton center based on perturbation theory for dark solitons. 



Introduction. Atomic Bose-Einstein condensates 
(BECs) constitute ideal systems for studying nonlinear 
macroscopic excitations in quantum systems flj. Excita- 
tions in the form of dark solitons and vortices are known 
to arise spontaneously upon crossing the phase transi- 
tion [1, 0] j a feature also studied in high-energy Q and 
condensed-matter [5( systems, in dynamical processes [jal 
and through controlled engineering 0, H, H, HH, Il2j |. 
In the latter category dark solitons are imprinted in a 
controlled manner after the gas has equilibrated @, [1, H, 
Hol . [TlT | . Although thermal effects revealed rapid soliton 
decay near the condensate edge @, HH, recent experi- 
ments at reduced temperatures (T <C 0.5T C ) [ij 
found the predicted [14] oscillatory pattern for the aver- 
aged soliton trajectories. 

To date, finite temperature dynamics of dark soli- 
tons have been investigated with phenomenological [l5| . 
quasiparticle scatteri ng J l6 fl . and generalised mean field 
[la ] models; see also [13, LL8| for quantum effects in var- 
ious background potentials. The former predict oscil- 
lations with increasing amplitude ('anti-damping' [01]), 
and appear to reproduce the average soliton trajectories 
to varying degrees of accuracy, however fail to account for 
the random nature of the experiments. In particular , ex - 
periments showed variations from shot to shot [^.[Tol.fll|. 
with single experimental realisations revealing the exis- 
tence of dark solitons for times much longer than those 
for which a reproducible (or average) pattern can begen- 
erated, an effect attributed to 'preparation errors' [9J. 

In this Letter we show that a spread in the trajectories 
of dark solitons prepared in the same manner could also 
arise due to the critical dependence of individual solitons 
on local phase/density fluctuations. Modeling the soli- 
ton dynamics by the Stochastic Gross-Pitaevskii Equa- 
tion (SGPE) [HUS] enables us to: (i) obtain an ab initio 
calculation of the spread of individual soliton trajecto- 
ries (Fig. 1, top); (ii) demonstrate that although aver- 
aging over different trajectories generates a well-defined 
pattern, this is restricted to times much less than the 
longest observed trajectories (Fig. 1, bottom), consistent 
with experimental findings [1, [lCj, El ; (hi) show that re- 



sults based on stochastic trajectory averaging can be well 
captured by the dissipative GPE (DGPE) [ll [U H [H , 
with an ab initio obtained damping coefficient; (iv) de- 
rive an analytical equation for the soliton center which 
captures such average dynamics very well at low temper- 
atures. 

Stochastic Dynamics: The SGPE [19, 20] describes the 
condensate and lowest excitations in a unified manner, 
including both density and phase fluctuations, with ir- 
reversibility and damping arising from the coupling of 
such modes to a thermal particle reservoir. Assuming a 
'classical' approximation for the mode occupations and a 
thermal cloud close to equilibrium, the SGPE reads [l|| 



ihdtip = (1 — vy) 



^-d 2 z +V(z)+g\^-n 



where g = 2ahio± is the effective ID coupling constant 
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FIG. 1: (Color online) Top: Normalised histograms of soli- 
ton decay times (main) and initial soliton depth, n so i, scaled 
to the average peak density (n(0)) (inset) (based on 200 re- 
alisations). Bottom: Individual stochastic trajectories from 
marked histogram bins (for as long as they are numerically 
tractable), 10-realisation trajectory average (black circles) 
and DGPE trajectory (green, dash-dotted). (Parameters: 

175nK, uj z = 2tt x 10Hz, 



N w 20 87 Rb atoms, T 
lo± = 2tt x 2500Hz, \v\ = 0.25c 
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(a is the scattering length, u>± 3> co z the transverse har- 
monic confinement), and V(z) = (l/2)mu> 2 z 2 the axial 
confining potential. 7 = i/3HI] K (z, t)/4 represents the 
a& initio determined dissipation arising due to the cou- 
pling to the thermal cloud (f3 — 1/k^T). £ K (z,i) is 
the Keldysh self-energy due to incoherent collisions be- 
tween condensate and non-condensate atoms and rj is a 
noise term with gaussian correlations (jf(z, t)rj(z' , t')) — 
2hk B Tf(z, t)S(z - z')S(t - 1') (see also 0, for further 
details and applications to condensate properties). 

Soliton experiments are modelled by first letting the 
system equilibrate at a given temperature and then in- 
troducing a dark soliton of specified velocity v in the trap 
center by multiplying ip by ip so \ = £tanh(£,z/£) + i(v/c), 
where C = yl — (v/c) 2 (£: healing length, c: speed of 
sound). Although the soliton generation is identical in 
all realisations (specified by v/c), fluctuations inherent 
in the atomic medium lead to a large variation in the im- 
printed soliton: The soliton speed v/c = \J\ — n ( i/n = 
008(5/2) is closely related to the depth of the density min- 
imum {rid) and the phase slip S across it. As a result, 
fluctuations in the background density upon generation 
should modify its depth, whereas the speed should also 
be affected by fluctuations in the condensate phase. 

The combination of these two factors leads to a slightly 
asymmetric spread in the initial soliton depth (Fig. 1, 
top inset), also interpreted as a stochastic change in 
the initial soliton speed (rs 30% for Fig. f). Moreover, 
the ensuing trajectory is further modified by the local 
phase/density fluctuations during the SGPE evolution. 

Soliton experiments are typically conducted in highly 
elongat ed g eometries, in order to avoid dynamical insta- 
bilities [23| . Phase fluctuations in such geometries set 
in at a characteristic temperature 7^ [2q ]. which can be 
much lower than the corresponding 'critical' temperature 
T c [13] • Although recent experiments were con- 

ducted in the regime T <C T$, T c , where both density and 
phase fluctuations are largely suppressed, soliton oscilla- 
tions can still be observed in the presence of phase fluctu- 
ations (T > 70 ), provided T <C T c . To amplify the differ- 
ences between individual trajectories, we thus choose re- 
alistic experimental parameters (N s» 20000 87 Rb atoms, 
uj z = 2t: x fOHz, lj±_ — 2h0oj z ) corresponding to this in- 
termediate regime 7^ <C T <C T c . This gives a phase 
coherence length w (O.f — 0.25)i? (R: Thomas-Fermi 
radius), with solitons allowed by ^> £. We focus on a 
relatively deep soliton which is more prone to this effect; 
we also anticipate phase imprinting to further enhance 
differences in trajectories due to the effect of fluctuations 
during the initial state preparation [l8| . 

Typical trajectories are shown in Fig. 1 (bottom) up to 
the point where the soliton can be numerically identified 
over the fluctuating background, which sets a decay time 
for each realisation. We find an asymmetric distribution 
of decay times, with some very long-lived trajectories. 
The spread in the decay times can be best visualised 
via characteristic trajectories from different histogram 
bins (labelled (a)-(c)). Despite their apparent differ- 
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FIG. 2: (Color online) Mean soliton decay times as a function 
of temperature obtained from the SGPE (red circles), or sin- 
gle DGPE realisations with ab initio determined 7(2) (green 
diamonds), or averaged 7 over [— R/2,R/2] (blue squares). 
Grey band indicates SGPE values within one standard de- 
viation of the mean decay time (for 200 runs); DGPE data 
shows time for soliton to decay to a depth comparable to 
the average background density fluctuations, with error bars 
corresponding to depths within the standard deviation of 
the fluctuations. Dotted horizontal line indicates time for 
one oscillation. Inset: 7(2) (solid) for T = 150 nK (bot- 
tom) and 300 nK with 7 (horizontal) for T=300nK; vertical 
lines show R — 281 z (parameters as in Fig. 1; characteris- 
tic temperatures in ID: = N(hui z ) 2 /ksfJ, « 25 nK [26| ]. 
T c = NhLo z /k B ln(2N) » 900 nK Q). 



ences, averaging over a sufficient number of trajectories 
(typically > 10) washes out such sensitivity, generating 
an antidamped oscillatory pattern, with a temperature- 
dependent shift in both amplitude and phase (black cir- 
cles). The average trajectory is only defined up to the 
earliest decay time within the set of trajectories con- 
sidered (here 27ui~ 1 ), in analogy to the experimentally 
reproducible soliton dynamics being restricted to much 
shorter times than those of individual long-lived trajec- 
tories d, [l(| [ll| . The average trajectory is practically in- 
distinguishable from an individual trajectory taken from 
the mean decay time bin (solid red), enabling us to in- 
fer the subsequent average soliton evolution from a single 
trajectory with a decay time close to the mean. 

Fig. 2 shows the dependence of the soliton decay time 
on temperature (red circles) in the intermediate temper- 
ature range of noticeable anti-damping: At higher T the 
soliton is lost to the fluctuating background, prior to ex- 
ecuting one full oscillation, thus leading to a decrease in 
the width of the decay time histogram, and to smaller 
error bars in the mean decay time; although our model 
predicts very little damping for T < 100 nK w 10%j;, 
consistent with recent pure condensate experiments [111 ] , 
our results may overestimate the actual lifetimes, due to 
the neglected role of collisions in the thermal cloud [l3| ■ 

The distribution of imprinted solitons and decay times 
is the main numerical result of this paper. Nonetheless, 
dynamics consistent with the average stochastic results 
can also be obtained by a simpler model discussed below. 

The DGPE and its comparison with the SGPE. A dis- 
sipative mean-field equation similar in form to Eq. |T]), 



3 



but without a noise term, was first introduced in a phe- 
nomenological manner by Pitaevskii [21]; in the BEC con- 
text this was applied to damping of excitations [23, vor- 
tex lattice growth [23l . l28j and dark soliton decay [l5| . A 
numerical advantage of the DGPE (which also restricts 
its predictive ability) is that only a single realisation is 
required under the assumption that trajectory-averaged 
properties should only depend on the dissipation; this is 
shown in Fig. 2 (inset): the self-consistent inclusion [29| 
of the mean field potential 2g(|?/>| 2 ) in the expression for 7 
generates a relatively flat profile around the trap center, 
with peaks at the condensate edges where the thermal 
cloud density is greatest. Since in the relevant dark soli- 
ton studies, the soliton spends most of its time well within 
the condensate, we can extract an averaged dissipation 
7, over a spatially-restricted region, e.g. 7 = J j(z)dz/ R, 
within [-R/2, R/2]. A simple analytical formula in the 
literature predicts 7(0) = a(ma 2 kBT /irh 2 ), with a « 3 
[HI . We find the spatially averaged rate 7 reveals a more 
pronounced scaling with temperature, though a reason- 
able first estimate can be obtained in the examined tem- 
perature range using this formula with 1/2 < a < 4. 

At low temperatures, the DGPE soliton oscillations are 
practically indistinguishable from the SGPE ones (Fig. 
1, bottom). A systematic comparison can be done by 
quantitatively comparing soliton decay times (Fig. 2): in 
the DGPE, these are identified by the time taken for the 
soliton to decay to a depth comparable to the background 
density fluctuations (as predicted here by a single SGPE 
run, or, in general, measured experimentally). We find 
very good agreement for both 7(2) and 7, within the error 
bars (grey bands), with a smaller relative error at lower 
temperatures. Since the DGPE reproduces the averaged 
results well in this regime (see also Fig. [3]below), we now 
provide an analytical solution for the soliton evolution. 

Analytical Results. Upon dropping the position depen- 
dence of 7(2) and further introducing the transformation 
£->(! + 7 2 )i, the ID DGPE takes the form: 



(i - j)d t ip 



-dl + V(z) + 



(2) 



where the density , length, time and energy are re- 
spectively measured in units of 2a, = yf h/muj^, cjJ 1 
and Hloj_, and V(z) = (l/2)f2 2 z 2 , with £1 = lu z /lu±_ < 1. 
We seek a solution of Eq. |2| in the form i/j(z, t) = 
^pb(z,t)e~ te ^v(z,t), where ipb(z,t) and 6{t) denote the 
background amplitude and phase respectively, while the 
dark soliton v(z,t) is governed by 



id t v + -d 2 z v - iPl{\v\ 2 - l)v 



d z j>i 



d z v + jd t v. (3) 



We assume that the condensate dynamics involves a 
fast scale of relaxation of the background to the ground 
state (justified a posteriori) and that the dark soliton 
subsequently evolves on top of the relaxed ground state. 
In the Thomas-Fermi limit, ip 2 sa fj,— V(z), and rescaling 
t — ► fit, z — > yfjiz, we obtain from Eq. a perturbed 



nonlinear Schrodinger (NLS) equation: 

id t v + ~d 2 z v-{\v\ 2 -l)v = P(v), (4) 
where P(v) stands for the total perturbation, namely, 



1 

2772 



dV 



2(1- Id 2 ) Vv + —d z v + 2-yud t v 
dz 



(5) 



and all terms in P are assumed to be of the same or- 
der (7 ~ fi). We now apply the perturbation theory for 
matter- wave dark solitons [30fl : starting from the dark 
soliton solution of the unperturbed system, we seek a so- 
lution in the form v(z,t) = cos (p(t) tanh?7 + i sin (p(t), 
where 77 = cos<p(i) [z — zq(£)], and ip(t) and zo(t) are the 
slowly- varying phase (\<p\ < tt/2) and center of the soli- 
ton. The resulting perturbation-induced evolution equa- 
tions for ip and zq, namely dip/dt = — (1/2) cos ipdV/dz + 
(2/3)7/iCOS</?sin<y9, and dzo/dt = sirup, lead to the fol- 
lowing equation of motion for the soliton center, 



d 2 z 
dt 2 



2 dzQ 
3 7 ^ 



Q V 
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dz 
dt 



(6) 



The nonlinear Eq. ([6]) can be integrated directly to yield 
the soliton trajectory: Fig. 3 shows very good agreement 
between the prediction of Eq. (6) (red) and the full DGPE 
(black) based on the spatially integrated 7, which are also 
consistent with the SGPE predictions with 7(2). 

In the case of a nearly black soliton (for dzo/dt suffi- 
ciently small), Eq. ^ is reduced to the linearized equa- 
tion d 2 z /dt 2 ~ (2/3)7Ai(dz M) + (n/V2) 2 z Q = 0. This 
includes the temperature-induced anti-damping term oc 
—jdzo/dt, and is reminiscent of the equation of mo- 
tion derived by means of a kinetic theory approach [l6| . 
For T = (7 = 0) the linearized equation recovers 
the constant amplitude oscillation of frequency £1/^/2 
0, [H. For T / (7 ^ 0), the solutions of the lin- 
earized Eq. ^ are Zo(t) oc exp(si^i), where s± t 2 = 
7At/3 ± \fK(ti/\/2) are the roots of the resulting char- 
acteristic equation. The discriminant A = (^i/^ cr ) 2 — 1 
(with 7cr = (3/n)(fl/V2) = 0.053 in our units) determine 
the type of motion: soliton trajectories are classified into 
sub-critical weak anti-damping (A < 0, 7 < 7 cr ), critical 
(A = 0, 7 = 7cr)j an d super- critical strong anti-damping 
(A > 0, 7 > j cr ) cases. Assuming an initial soliton loca- 
tion ^o(0) = and velocity io(0), the sub-critical soliton 
trajectory reads: 



z {t) 



io(0) 



aJ 



n / 7 2 

cos(w i), L0 o — — a 1 —, (7) 



V2 ] 



1c 



indicating an exponential increase in its maximum am- 
plitude (Fig. 3 top, dashed green line), whose magnitude 
depends on both temperature and chemical potential; the 
oscillation frequency lu is also shifted from its T = 
value [131 ]. Corresponding trajectories in the critical and 
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FIG. 3: (Color online) Top: Soliton DGPE trajectories 
(black) vs. analytical results (nonlinear: solid red; linear: 
dashed green) with 7 = 7 = 0.00014, and stochastic ones 
with 7(2) (stars: 10-trajectory average, circles: single mean- 
bin trajectory). Bottom: Dependence of the real part (left; 
instability growth rate) and imaginary part (right; oscillation 
frequency) of the unstable eigenmode of the excitation spec- 
trum on 7: analytical results (solid red) vs. DGPE numerics 
(blue circles) (Parameters as in Fig. 1, with T = 150 nK.) 



super-critical cases read: zo(t) — io(0)t exp(7/z£/3) and 
Zo(t) = [i (0)/(si - s 2 ))[exp(s 1 t) - exp(s 2 £)]. 

The above results are also supported by a linear sta- 
bility analysis around the stationary dark soliton, tpds- 
This waveform makes the right hand side of Eq. ^ van- 
ish and is, thus, an exact solution of the T ^ problem. 
As rigorously proven [3l|, the anomalous (or negative 



Krein signature) mode of the dark soliton leads to an in- 
stability, upon dissipative perturbations. In particular, 
the relevant mode of the linearization around the soli- 
ton (solution of the eigenvalue problem arising from ip — 
ipds + e(exp(\t)a(x) + cxp(X*t)b i '(x)) for the eigenvalue- 
eigenvector pair {A, (a, b)}) acquires Re(A) > for 7 > 0. 
Fig. [3] (bottom) demonstrates an excellent agreement be- 
tween the analytical prediction for the relevant eigenvalue 
and the numerical result for the excitation spectrum of 
the DGPE. 

Discussion: A full description of the rich experimental 
features observed in dark soliton experiments requires a 
stochastic model incorporating density and phase fluctu- 
ations, in order to model experimantally relevant shot- 
to-shot variations. The stochastic Gross-Pitaevskii equa- 
tion was shown to capture these features well, leading to 
specific predictions for the spread of soliton decay times 
with different realisations of the dynamical noise, in close 
analogy to different experimental realisations. Nonethe- 
less, even within the phase-fluctuating regime, mean soli- 
ton trajectories/decay times are captured reasonably by 
the simpler dissipative Gross-Pitaevskii equation (with 
additional experimental or theoretical input required to 
obtain the dissipation term). A fully analytical solution 
of the dark soliton motion in excellent agreement with 
the dissipative Gross-Pitaevskii equation applied to the 
finite temperature was given, paving the way for future 
analytical studies of other macroscopic excitations, such 
as vortices, in atomic condensates. 
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